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ABSTRACT 

We investigate a purely stellar dynamical solution to the Final Parsec Problem. Calactic nuclei 
resulting from major mergers are not spherical, but show some degree of triaxiality. With iV-body 
simulations, we show that massive black hole binaries (MBHB) hosted by them will continuously 
interact with stars on centrophilic orbits and will thus inspiral — in much less than a Hubble time — 
down to separations at which gravitational wave (GW) emission is strong enough to drive them to 
coalescence. Such coalescences will be important sources of GWs for future space-borne detectors such 
as the Laser Interferometer Space Antenna (LISA). Based on our results, we expect that LISA will 
see between ~ 10 to ~ few x 10^ such events every year, depending on the particular MBH seed model 
as obtained in recent studies of merger trees of galaxy and MBH co-evolution. Orbital eccentricities 
in the LISA band will be clearly distinguishable from zero with e > 0.001 — 0.01. 

Subject headings: black hole physics — galaxies: nuclei — stellar dynamics — gravitational waves 



1. INTRODUCTION 

Massive black hole binaries (MBHBs) are one of the 
most interesting sources of gravitational waves (GWs) 
for future space-borne detectors such as the Laser Inter- 
ferometer Space Antenna (LISA). They are expected to 
coalesce under the strong emission of GWs, after stellar- 
and/or gas-dynamical processes bring them to separa- 
tions small enough {acw ~ 10~^ pc) that GW emission 
is effi cient in making them coalesce in less than a Hubble 



time ( Milosavljevic fc Merritt|2003 Armitage & Natara- 



jan 20U5 T It is still an open problem whether MBMB 



coalescences are generic and prompt, or whether long- 
lived binaries are the norm. 

The paradigm for MBH binary evolution, after a 
merger of gas-poor galaxies, consists of three distinct 



phases (fBegelman, Blandford & Rees 1980). First, the 
two MBHs sink towards the center due to the dynamical 
friction exerted by the stars. This process continues after 
they form a bound pair at a semimajor axis separation 
a r/i, where r^ is the binary's influence radius defined 
to be the radius which encloses twice the mass of the bi- 
nary in stars. It stops when the binary reaches the hard 
binary separation a ^ Oh ( (Quinlan^^l996; ,Yu^,2002^ 
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where is the binary's reduced mass, cr is the local 
ID velocity dispersion, q = M,^2A^»,i is the binary's 
mass ratio. Secondly, for a < a^, as dynamical friction 
becomes inefficient in further driving the inspiral, it is in- 
stead the slingshot ejection of stars, following three-body 
scattering with the binary, that dominates. Thirdly, the 
binary eventually reaches a separation acw at which the 
loss of orbital energy to GW emission drives the final 
coalescence. The transition from the first to the second 
phase is prompt provided that the mass r atio of the rem- 
nants is not too small q = M^/M-] > 0.1 (Colpi & Dotti 
Callegari et al.|[20TT ). In contrast, the subsequent 
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transition from the second to the third phase could con 
stitute a bottleneck for the binary evolution towards final 
coalescence. This is the so-called Final Parsec Problem. 

In quasi-steady spherical stellar environments, the bi- 
nary's hardening rate s{t) = d/dt{\/a) slows down sig- 
nifica ntly once it reaches separations a few times below 
^ a u (iQuinlan &: Hernquist 1997[ Milosavljevic & Mer- 
|ritt| [2003; Bcrczik et al.| |2003| ^ In these spherical and 
gas-poor nuclei, two-body relaxation is the only mecha- 
nism for populating the binary's loss cone[^ but being a 
slow diffusive process, it is only in low-luminosity galax- 
ies harboring MBHs of mass M, < few x 10^ that 
central relaxation times are short enough to driv e the bi- 
nary to coa lescence in less than a Hubble time (Merritt 



et ah 20071. 



But spherical models are a worst case scenario — and 
not a very realistic one at that! Merger remnants will 
generally be irregular with some degree of triaxiality 
and, even if triaxiality would only be a rather mild and 
transient phenomeno n, it may suffice to bri ng the bi- 
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(2006) and 



nary d own t o agw (Merritt fc Poon 2004). 



ySSm) studied tri 
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Berentzen et aE 
rotating models ot galactic nuclei using TV-body 
simulations — including the full post-Newtonian correc- 



^ The loss cone is the region of phase space corresponding, 
roughly speaking, to orbits that cross the binary, i.e. wi th angular 



momentum J < Ji^ 
,fc Shapiro|1977| . 



\/GMi2fai,i„, where / = 0(1) ( Lightman| 
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tions to the MBHB. They have shown that MBHBs in 
such models do indeed coalesce in much less than a Hub- 
ble time. The next logical step is to study mergers of 
galactic nuclei to investigate if the latter results still hold 
true under more realistic models and initial conditions. 

In this Letter, we use N-body simulations to show 
that: (1) in merging nuclei, the hardening rate is N- 
independent — allowing the extrapolation of TV-body re- 
sults to real galaxies; (2) the triaxiality depends on the 
orbital parameters of the progenitor galaxies: prolate 
shapes occur when the merger is almost radial, while an 
oblate morphology is the result of a less radial merger; 
(3) MBHs become bound with high eccentricities (up 
to e ~ 0.95); (4) the eccentricity tends, on average, 
to incr ease in good agreement — often quantitative — with 
[Quinlan (1996) predictions; (5) high eccentricities assist 
the MBHB into promptly coalescencing in much less than 
a Hubble time; (6) eccentricities in the LISA band are 
likely to be distinguishable from zero (e > 0.001 — 0.01) 
even though GW circularizes the orbits, and will also be 
quite large (0.4 < e < 0.8) in the Pulsar Timing Array 
(PTA) band. 

2. MODELS AND INITIAL CONDITIONS 

We have performed two sets of iV-body experiments. 
In both sets, galactic nuclei ar e represented by spheri- 
cally symrn etric Dehnen models ( Dehnen|1993 Tremaine 
et al.|1994| ) . These models have a central power law den- 
sity profile, p{r) = (3 — j)MT/4:TTr'^ {rt + f)^~^, with 
logarithmic slope 7 and a break radius r^, which are both 
set equal to one. The total mass of each nucleus is set 
Mt = 1, and we adopt units where G — 1. The total 
mass of the binary MBH M12 — M,^i + M,_2, and we 
take q = M,, 2/Af,, 1 = 1. We study unequal mass MBH 
coalescence s in parallel papers ( Berczik et al. 12011 Preto 
erar||20llj ). 



I'he set (A) consists of a single spherical nucleus where 
two MBHs are placed symmetrically about the center, on 
an unbound orbit, with initial separation Aro = 2, initial 
angular momentum L/ Lc = 0.5, where Lc is the angular 
momentum of the local circular orbit. The set (B) con- 
sists on the equal-mass merger of two initially bound — 
but well-separated — spherical nuclei, each of which has a 
single MBH at the center with zero initial velocity with 
respect to its nucleus. For B, the initial separation Aro 
refers to both nuclei taken as if they were point masses 
located at each center of mass. The half-mass radius of 
each nucleus is ri/2 ~ 2.41; accordingly, and in order 
to have an initial configuration with two well separated 
nuclei, while minimizing the computing time, we set the 
initial separation equal to 20. For the initial orbital an- 
gular momentum of the binary nuclei, we have taken two 
values L/Lc ~ 0.14 and 0.6 given the nearly-parabolic 
encounters typical of major galaxy mergers seen in cos- 
mological simulations ( Khochfar fc Burkert| 2006 ). Dur- 
ing the first pericenter passages, the MBH separations 
are Atbh 0.2 ~ O.lri/2 and Atbh ^ 2.2 ^ ri/2, 
respectively. Table 1 lists the runs and adopted param- 
eters. 

We have performed the A^-body simulations using the 
parallel yj-GPU code. This is a yet unpublishe d vari- 
ant of the parallel direct N-body code (y9-GRAPE ( Harfst 



et al. 200 7|, which uses GPU accelerator cards on parallel 
clusters (iBerczik et al.||2011|). It includes a fourth-order 
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TABLE 1 

af-body integrations. 1^* column: mass of the mbh binary; 
Other columns: particle number N; First two lines: 

simulations of spherical nuclei; last four lines: 
Simulations of merging nuclei; A = L/Lc measures the 

INITIAL ORBITAL ANGULAR MOMENTUM OF THE MBH BINARY 
(Asph = 0.5 FOR SPHERICAL NUCLEI), OR OTHERWISE IT MEASURES 

THE INITIAL ORBITAL ANGULAR MOMENTUM OF THE MERGING 
NUCLEI (Al = 0.14 FOR NEAR-RADIAL MERGER AND A2 = 0.6 FOR 
LESS RADIAL MERGER). AlL NUCLEI HAVE 7=1; ALL BINARIES 
HAVE EQUAL MASS q = M, ^2/^,^1 = 1. 



Hermite integration scheme, with b lock time steps, anal- 
ogous to NBODYl ( |Aarseth||2003| . 

The code does not include regularization of close en- 
counters, and softening of the gravitational interaction 
is adopted instead. The softening length has to be cho- 
sen small enough that it reproduces the refilling of the 
binary's loss cone by two-body relaxation. After some 
testing, we adopt a softening length e = 10""^ in model 
units. We set the time step parameter ( Aarsethj 2003 1 
to 77* — 0.01 for the field stars and tibh = O.UOl tor the 
BHs. Furthermore, we force the MBHs to be advanced 



synchronously at all times with the smallest step. With 
the parallelized version of the (^-GPU code, one can study 
models with very large nu mber N of parti cles and the re- 
sults agree with NB0DY4 ( Aarseth|2003 ) as far as single 
stars and distant encounters are concerned. For the high 
velocity dispersions present in nuclei with a MBH, the 
effect of close encounters between field stars is negligi- 
ble for the bulk evolution of the stellar system ( |Preto fc] 
Amaro-Seoane||2010[ ). 



3. MBH EVOLUTION IN SPHERICAL VERSUS IN 
MERGING NUCLEI 

The stars that drive the orbital decay of a hard MBHB 
are those that enter the loss cone orbits. The MBHB's 
hardening rate is thus determined by the product of the 
flux of stars entering the loss cone with the average ki- 
netic energy they receive when ejected — at the expense 
of the MBHB's orbital energy — through the slingshot 
mechanism. Denoting by Tic{E,t) the time-dependent 
flux into the loss cone and by {AE{E)) the mean kinetic 
energy imparted to stars which are sc attered off ' by the 
binary, the hardening rate is given by ( |Yul[2002l ) 



d_ 
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dE{AE{E))Ti,{E,t), (2) 



where E = GMi2/r + <i>*(r) - 1/2 i;^, and <i>*(r) is the 
gravitational potential due to the stars. The mean ki- 
netic energy {AE{E)) is given by 



{AE{E)) ^ {C) 



G^r 



(3) 



where (C) ~ 1.25 is a dimensionless quantity whi ch was 
measured from three-body scattering experiments dQuin- 



|lan^,1996J . Therefore, the hardening rate s{t) can be 
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Fig. 1. — Binary hardening. Upper panel: in a spherical nucleus, 
s{t) decreases with A^. Middle panel: in a merging nucleus, s{t) is 
N-independent. Lower panel: hardening rates as a function of N 
for different M12. Being much smaller, (s) of the M12 = 0.1 binary 
has been multiplied by 100 to better fit in the plot. Labels 's' for 
spherical and 'm' for merger. 



rewritten as 

s{t) - ^ 

dt \ a 



2m,{C) 
Ml 2 a 



+ 00 



dETic{E,t). (4) 



The time evolution of the flux J^ic{E,t) depends on 
the symmetries of the gravitational potential — and on 
the orbit families it supports — of the nuclei in ques- 
tion. In principle, J-ic{E,t) in the spherical case can 
be obtained from Fokker-Planck calculations that take 



into account the diffusion of stars in phase spac e (Mer- 
ritt et al. I [20071 [Preto &: Amaro-Seoaiie| |2010[ ). Here 
we derive simple scaling relations which are useful in 
interpreting the N-body results. For each energy E, 
Tic{E,t) oc n{E,t)/Trix{E,t) where n{E,t) is the num- 
ber of stars of energy E per unit energy and Trix{E,t) 
is the local two- bodv rel axation time; the latter scales as 
Trix OC I pnii, ( Spitzerjl987 ). The flu x of stars into the 
loss cone is expected to peak around Th ( Perets & Alexan- 
dCT|[2008| , so we evaluate these quantities there. Hence, 

f Afi2)/rft ~ 3GMi2/r^ cx M\!^^— 

follo ws from the Af, — a relation 

2005|. Then, a\ cx Aff^** obtains. 



al ~ G(Af(< r^) 
where cx Ai", 



1/2 



(Ferrarese & Ford 



On the other hand, tor a fixed galaxy mass, we have 

cx 1/iV and therefore T^ix cx m\!^N j p. Since in our 
N-body models, p{r) and n{E^ t) are unchanged and only 
a changes as A/12 is varied, we find that the hardening 



Fig. 2. — Triaxiality T and mass flattening e = (a — c)/a of merg- 
ing nuclei. Shown are merging binaries of total mass M12 = 0.005 
with L/Lc = 0.6 (upper panel), M12 = 0.01 and 0.02 with 
L/Lc = 0.14 (almost radial mergers, in the middle and lower pan- 
els). T and e are measured in five mass shells between r = and 
r = 2.5, each of width Ar = 0.5. Triaxiality decreases over time, 
the faster the heavier the binary is. Mass flattening is constant. 

rate scales with M12 and as s cx Afj"//*A^-i. The case 
of a triaxial nucleus is different: J for each star is not 
conserved, thus stars m ay precess into the los s cone on 



a time scale Tpr < T^-ix ( |Merritt fc Poon||2004[ ); and Tj 



will depend only on the global gravitational potential of 
the galaxy. In this case, the mass flux into the loss cone 
m^J-iciEji) cx m^n{E,t)/Tpr{E,t), and also s(i), will be 
independent of the number N of stars. 

In Figure [1] we see that s{t) is N-dependent in a spher- 
ical nucleus, while it is A'^-independent in the merging 
one. In the former case, s{t) cx N~", with a = 0.45 and 
0.75 for binaries of M12 = 0.005 and 0.1. These resuhs 
can be interpreted as follows. In the empty loss cone 
limit {a = 1), the stars repopulate the loss cone at a rate 
cx T^i^ much lower than the that with which they are 
ejected by the binary, which is cx 

Tdyi- In the fuU loss 
cone limit {a — 0), stars enter the loss cone at a rate 
which is similar to the rate at which they are ejected by 
the binar y. A measure for the loss c one refilling rate is 
given by ( jLightman fc Shapiro||1977| 



qiE) 



6J_ 

Jlc 



(5) 



where SJ is the mean change in J, per orbital period, 
of a star on a low- J. In the limit when q{E) <C 1, the 
loss cone is said to be empty; while q{E) >• 1 in the 
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Fig. 3. — Long term eccentricity evolution. Red and green 
lines represent NB and semi-analytic evolution without radiation 
reaction. Blue and magenta lines correspond to semi-analytic 
solution, including radiation reaction, for M12 = IO^Mq and 
M12 = 10* Mq, respectively. 

full loss cone limit. For a given nucleus, and for r > r^, 
we expect 5 J to be independent of M12. As a result, 
the weaker dependence of (s) on N for lighter binaries, 

placed in a spherical nucleus, follows from q cx M^2^'^') 
at the same E, q{E) of the M12 = 0.005 binary is ^ 4.5 
larger than that of the M12 = 0.1 one. We would need to 
use a larger N for the lighter binaries, (s) cx A^^'' '*^, to 
enter deep into the empty loss cone limit (s) cx N^^; the 
heavier binary, (s) cx iV~"-^^, almost reaches this limit. 
The dependence of (s) on M12 is more straightforward 
to interpret. For the spherical case, the lighter binary 
is expected to harden at a rate ^ 20'^/'' higher than the 
heavier, which is indeed the case. In the merger case, 
m^Fic{E) is A^-independent and therefore (s) cx M^^ ■ 
Since the mass ratio between the binaries is 2, (s) also 
differs by a corresponding factor of two. 

Following [Merritt fc Poon| ([2Cl04| , we measure the tri- 
axiality of the nucleus with T = [a? — b"^) / {a? ~ c^) . 
Figure [2] depicts the evolution of T and of the flatten- 
ing e = 1 — c/a for several mass shells of merging nuclei. 
The value of T of each remnant, immediately after the 
merger, depends on the initial L/ L^.. In the case of a near 
radial merger, L/L^ = 0.14, the remnant is prolate and 
evolves over time towards an oblate spheroidal shape; 
for L/Lc = 0.6 the remnant is an oblate spheroid from 

* Models with T = 0.25 and 0.75 correspond to moderately 
oblate and prolate shapes, respectively. 
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Fig. 4. — Upper panel: Range of coalescence time for binaries 
with Mi2 = 1O''M0 and 10^ Mq. Middle panel: distribution of 
eccentricities, for M12 = 10^ Mq, at ai,i„ = lOORschw Lower 
panel: distribution of eccentricities, for M12 = IO^Mq, when 
/orb = Vgw = 10-«Hz for PTAs. 

the very beginning. The triaxiality decreases over time, 
and the rate at which it changes is faster the larger M12 
is. The triaxiality remains significant, in the inner mass 
shells, until the binary reaches the relativistic phase in all 
models with the smallest (and more realistic) values of 
Afi2, and also in most of the other cases. The flattening 
e ^ 0.2 is constant throughout in all cases, so the asymp- 
totic shape of the merger is that of an oblate spheroid. 
We conclude that the rather mild triaxiality created dur- 
ing the merger supports a family of centrophilic orbits 
that keep the loss cone full {a = 0) at all times until the 
binary reaches relativistic separations ^ o-gw- 

4. ECCENTRICITY EVOLUTION AND TIME SCALES FOR 
COALESCENCE 

The hardening rate, due to the slingshot ejection of 
stars, is found to be in our A^-body simulations es- 
sentially independent of the binary's orbital elements, 
while that due to GWs is strongly dependent on them: 
d/dt( l/a)Gw ~ Id/a^law oc a'^il - e^)-^/^ ( [Peters 
19641. As a result, the time a binary takes to coalesce 
depends strongly on its eccentricity. In paper I, we did 
follow this evolution self-consistently with N-body simu- 
lations of rotating King models. Since such calculations 
are extremely CPU-intensive, we est imate the fu ll evolu- 
tion using a semi-analytic approach ( Quinlan|19 96). The 
advantage is that we can calibrate the average hardening 
rate (s) with our N-body simulations — which would re- 
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main a free parameter otherwise — , and thus make quan- 
titative predictions on both the coalescence times and the 
long term eccentricity evolution. 

The evolution of the MBHB orbital elements, including 
the effect due to orbital energy lost to GWs, is given by 



di 



d 
di 



d 
di 



GW 



de 
di 



de\ 

/ GW 



(6) 



The GW terms are given in Peters (1964). The eccen- 
tricity evolution, driven by the stars, is obtained from 
three-body scattering experiments 



^K{e) a {s), 



(7) 



where K(e) = e(l — e^)''°(fc i-|-fc2e) and the constants are 
taken from Quinlan (1996). In order to assess the qual- 



ity of the fits to the iV-body results, Figure [3] compares 
the A^-body evolution of the binary with that obtained 
from the semi-analytic model. We take as initial condi- 
tions for the integration of equations (ro| an instant of 
time in the early hard binary phase. Given the differ- 
ences between our A''-body models and the assumptions 
embodied by the semi-analytic description the agreement 
is quite remarkable. 

We then include the GW terms due to radiation reac- 
tion to compute the time it takes for the binary to coa- 
lesce. To scale our models to binaries with M12 = 10^ 
and M12 = 10® Mq, we adopt the most recent observa- 
tional values for the mass normalization of the Milky 
Way nucleus, M(< Ipc) = IO^Mq ( |Sch5del et aL||2009 ) 
and use the M, — a relation to extrapolate to dirterent 
MBH masses. The results are shown in the upper panel 
of Figure |4l We see that coalescence times range between 
Tcoai ~ 10^ yrs and ~ few x 10* yrs. These times are 
not longer that the mean time between successive major 
mergers. In contrast, for a spherical nucleus, coalescence 
times for the lower mass would become ^ few x Gyr, 
while binaries with > 10* M0 would stall (Preto et al. 



2011|. 

We also follow the long term evolution of the eccentric- 
ity. In the TV-body runs, the binaries become bound with 
high eccentricities (up to e ^ 0.95) on average — in agree 



ment with previous works (Berentzen et al. 2009 Preto 



et al.| 2009). Since LISA will be sensitive to the inspiral 



signal of IO^Mq binaries, it is important for data analy 
sis purposes to estimate whether they will ente r the band 
with non-ne gligible eccentricity (e > ^^Zt^ ( [Porter fc 



SesanapOlO ). The middle panel of Figure |4] displays the 
distribution of eccentricities at a = lOO-Rsc/iuQ — most 
binaries will not be fully circularized by then. We ex- 
pect therefore that the eccentricity in the LISA band will 
be non- negligible. Finally, the lower panel of Figure [4] 



depicts the eccentricity distribution at few = "^forb — 
10~*Hz for the PTA band. We see that eccentricities 
are quite high — peaking at e ~ 0.6. The results pre- 
sented here concerning the coalescence times and eccen- 
tricity growth corroborate recent three-body scattering 
studies — which had to t reat the aver age hardening rate 
(s) as a free parameter ( Sesana|[2010 ). 



5. SUMMARY 

With our results, we are moving closer towards a con- 
sistent solution to the Final Parsec Problem, and thus 
of providing a dynamical substantiation to the cosmo- 
logical scenario where prompt coalescences are assumed 
during the course of galaxy evolution (Sesana et al. 2007j 
Volonteri||2010 ). Our results suggest that the formation 
of eccentric binaries, followed by a quick orbital decay, 
could result from the expected development of global 
non-axisymmetries in galaxies after they merge. Our gas- 
poor merger models show only rather mild departures 
from axisymmetry and a small amount of rotation; we be- 
lieve that stronger departures from axisymmetry — to be 
expected from higher amount of rotation — , and the pres- 
ence of gas will only reinforce our conclusions. It seems, 
therefore, probable that prompt coalescences result from 
mergers of irregular galaxies expected to be common at 
high redshift. Based on our prompt MBHB coalescence 
results, we expect that LISA will see 10 — few x 10^ 
events per year depending on the MBH seed model 
(Sesana et al. 2009 Volonteri 2010). Moreover, even 
though GWs circularizes the MBHBs during the late rel- 
ativistic phase of inspiral, they are likely to have some 
residual (e > 0.001 — 0.01) eccenticity when entering the 
LISA band and a broad distribution (0.4 < e < 0.8) in 
the PTA band. 
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